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Abstract 

In this paper we derive control charts for the variance of a Gaussian process 
using the likelihood ratio approach, the generalized likelihood ratio approach, the 
sequential probability ratio method and a generalized sequential probability ratio 
procedure, the Shiryaev-Roberts procedure and a generalized Shiryaev-Roberts ap- 
proach. Recursive presentations for the calculation of the control statistics are 
given for autoregressive processes of order 1. In an extensive simulation study these 
schemes are compared with existing control charts for the variance. In order to 
asses the performance of the schemes both the average run length and the average 
^ delay are used. 
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o 1 Introduction 

In many applications we are faced with the problem to detect changes over time in an 
observed process. Because usually a change should be detected fast sequential methods 
t-h are more appropriate in such a situation. The most important tools for monitoring a 

process are control charts (cf. Stoumbos et al. (2000)). Control charts are successfully 
applied in engineering for a long time (e.g., Lawson and Kleinman (2005), Frisen (2007)). 
In the last 20 years many further applications in different areas have been studied like, 
e.g., in public health, economics, environmental sciences. In that context the underlying 
processes have a more complicate structure and are mostly modeled by a time series. 
Alwan and Roberts (1988) showed that control charts for independent variables cannot 
be directly applied to time series. They proposed to use residual charts, i.e. to transform 
the original observations such that the transformed observations are independent and to 
apply the well-known control charts to these residuals. Residual charts have been studied 
by several authors (e.g., Harris and Ross (1991), Wardell et al. (1994), Lu and Reynolds 
(1999)). Another possibility is to directly monitor the observed process. The behavior 
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of the Shewhart chart for time series was studied in Schmid (1995). An extension of the 
exponentially weighted moving average (EWMA) chart of Roberts (1959) to time series 
was proposed by Schmid (1997a). Cumulative sum (CUSUM) charts for time-dependent 
processes have been studied among others by Nikiforov (1975), Schmid (1997b), Frisen 
and Knoth (2012). An overview on control charts for time series is given in Knoth and 
Schmid (2004) and Okhrin and Schmid (2007). 

Most of the literature on that topic is dealing with the monitoring of the mean behavior 
of the observed process. Here we want to focus on the surveillance of the variance of a 
time series. We are interested to detect an increase in the variance. For instance, such 
a question is of great importance in economics where the variance is the most applied 
measure for the risk and an early detection of a change in the risk behavior of an asset 
is an important information for an analyst. The first EWMA control chart for time 
series was introduced by MacGregor and Harris (1993). Schipper and Schmid (2001) 
introduced several one-sided variance charts for stationary processes, however, their main 
focus was in the area of nonlinear time series. In many applications the control statistic 
for independent processes is used and the independent process is replaced by the time 
series. Thus the structure of the time series is not taken into account for the derivation 
of the control statistic. This is a disadvantage of this procedure. 

In this paper we derive control charts for the variance of a Gaussian process by making 
use of the likelihood ratio approach (LR), the sequential probability ratio test of Wald 
(SPRT), and the Shiryaev- Roberts (SR) procedure. For deriving these charts it is first 
assumed that the size of the change is known. Thus all obtained charts depend on a 
reference value which has to be suitably chosen in advance. This is sometimes a drawback 
in applications. We consider generalized control charts as well. They are obtained via 
the generalized LR, SPRT, and SR procedure. The great advantage of these schemes is 
that they do not depend on a reference value. It has to be emphasized that our results 
are quite general and cover all autoregressive moving average processes with a Gaussian 
white noise. 

In Section 2 the underlying model of the paper is introduced. It is explained how in 
our paper the target process and the observed process are related with each other. In 
Section 3 the CUSUM control chart for the variance in the independent case is briefly 
presented and a new CUSUM control scheme for Gaussian processes is derived over the 
LR approach. In an example we consider the special cases of autoregressive processes of 
order 1 and 2 and it is shown that in that case the control statistics can be calculated 
recursively. In Section 4 a CUSUM variance chart is derived by the SPRT and the result 
is a residual chart. In Section 5 the Shiryaev- Roberts method is used to get a control 
chart for the variance. In the Sections 6 to 8 generalized control schemes are derived. In 
Section 6 we use the generalized LR method, in Section 7 a generalization of the SPRT 
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approach, and in Section 8 a generalization of a modified version of the SR method is 
obtained. 

In an extensive simulation study these control schemes are compared with each other 
assuming that the underlying target process is an autoregressive process of order 1 (AR(1)) 
(Section 9). As a measure for the performance of a control scheme the average run length 
(ARL) and the average delay are taken. All charts are calibrated such that the in-control 
ARL is the same if no change is present. Our results show that except the SR chart all 
other schemes with a reference value have the smallest out-of-control ARL if the reference 
value is equal to the true value of the change. It turns out that the generalized SR chart 
has the smallest ARL if the change is small. It is even better than the charts with the 
optimal reference value. For medium and larger changes the LR chart and the SPRT chart 
provides better results provided that the reference value is not dramatically smaller than 
the true change. Except the generalized SPRT scheme for all charts the worst average 
delay is equal to the average run length. The limit of the average delay seems to be 
the smallest for the GSR chart. This scheme must be preferred for larger changes if the 
change arises at a later time point. 



The aim of statistical process control (SPC) consists in detecting structural deviations 
within a process over time. It is examined whether the present observations can be 
considered as realizations of a given target process {Y t }. The procedure is a sequential 
one. The observations (samples) are analyzed consecutively. It is desirable to detect a 
change as quickly as possible after its occurrence. Of course there are various types of 
changes which may influence the target process. In this paper we focus on the detection 
of an increase in the variance. Such a problem arises in practice very often. For instance, 
the variance is considered as a risk measure in economics and thus an increasing variance 
is a hint that the risk of an asset is getting larger. In engineering the variance reflects 
the quality of a production process and an increase is a bad sign since the production is 
getting worse. 

Let {Yt} be a (weakly) stationary process with mean fi and autocovariance function 
Cov(Y t , Y t+ h) = l{h). In what follows it is assumed that the relationship between the 
target process {Y t } and the observed process {X t } is given by 



for t G Z with A > 1 and tGNU {oo}. Thus a change in the scale appears at position 
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for 1 < t < t 



(1) 
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r if r < oo. The observed process is said to be out of control. Else, if r = oo, then 
{X t } is called to be in control. Here it is assumed that at a given time point exactly one 
observation is available. 

Note that the change in the scale does not influence the mean structure. It holds that 
E{X t ) — (j, as well in the in-control state as in the out-of-control state. Moreover, we get 
that 



Var{X t ) 



Cov{X u X t+h ) = { 



jvar(Y t ) for l<t<r 

[A 2 Var(Y t ) for t>r 

7h for t < min{r, r — h} 

A'jh for min{r, r — h} < t < max{r, r — h} ■ 

A 2 7 ft for t > max{r, r — h} 



Thus the observed process {X t } is not stationary in the out-of-control case. 

3 A Variance Chart based on the Likelihood Ratio 
Approach 

3.1 LR Approach applied to Independent Variables 

The CUSUM chart for the mean was introduced by Page (1954). The CUSUM scheme 
is a control chart with memory. At each time point the decision is based not only on 
the present observation, but also on previous observations. The CUSUM chart gives 
all former observations the same weight. Moreover, the control chart depends on an 
additional design parameter, the reference value. Up to now mainly CUSUM charts for 
the mean have been considered in literature. 

CUSUM charts can be derived by means of the log likelihood ratio approach (see 
Siegmund (1985, Ch. II. 6)). In that context it is demanded that A is a known quantity. 
For instance, assuming that the variables Y t are independent and normally distributed 
with expectation /i and variance 7 the log likelihood ratio of the present model is given 
by (see, e.g., Hawkins and Olwell (1998) and Schipper and Schmid (2001)) 

1-1/A 2 ~ 



for n > t with 



S n (A) = j2 i -^- J ^-nK(A) (2) 
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and 

We conclude that the process is out of control at time n > 1 if 

Si (A) = S n (A) - mm &(A) > c. (4) 

0<t<n 

Note that S+(A) can be calculated recursively as S^(A) = max{0, S'^_ 1 (A) + (X n — 
/ i ) 2 /7o — -^(A)} for n > 1 with Sq(A) = 0. This representation dramatically simplifies 
the practical calculation of the control design. 

Unfortunately A is unknown in practice and for that reason it is necessary to fix a 
value for A, say A*. In practice A* > 1 is interpreted as the change against which we 
want to be protected. Then we have to replace K(A) by K(A*). Thus we obtain 

S+(A*) = max{0X-i(A*) + ~ ^ - K(A*)},n > 1,S+(A*) = 0. (5) 

7o 

The process is concluded to be out of control at time n > 1 if S„(A*) > c. For each 
fixed value of A*, we determine the value of c such that the in-control ARL is equal to 
some predetermined quantity £. 

In the above derivation it is assumed that the underlying process is independent. In 
this paper, however, we are interested in target processes which follow a time series. 
Because of the complicated structure of the likelihood function of a time series practition- 
ers have used the above recursion (|5| for time series as well (e.g., Schipper and Schmid 
(2001)). Then {X t } stands for a time series, \x is equal to its mean and 70 is the in-control 
variance. It has to be noted that the control limit c must be chosen by taking the time 
series structure into account. The run length of this scheme is given by 

N LR , ud (c; A*) = un> e N : S+(A*) > c} (6) 

with S+ (A*) as in (g). 

While for independent observations several numerical methods are available to cal- 
culate the ARL (cf. Brook and Evans (1972), Knoth (2010)) the determination of this 
quantity turns out to be very difficult if the target process has a dependence structure. In 
that case no explicit formulas for the ARL are available and in practice the control limits 
are determined via simulations. 
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3.2 LR Approach applied to Gaussian Processes 

In this section {Y t } is assumed to be a Gaussian process with mean zero and covariance 
function k(i,j) = Cov(Xi, Xj). Assume that the covariance matrix [k(i, j)]ij=i,.., n is non- 
singular for each n > 1. Let Y t denote the best linear predictor of Y t based on Y t _i, ..,Y 1 . 
This quantity can be recursively calculated using the innovations algorithm (cf. Brockwell 
and Davis (1991, p. 172)) 



for t = 1 



The quantities 9tj can be determined recursively Note that Y t = X]*=i a tjYj with some 
coefficients a t j. Note that for a Gaussian process the best linear predictor is equal to the 
predictor obtained by minimizing the mean-square distance. 

Following Brockwell and Davis (1991, p. 255) the likelihood function of (Yi,..,Y n ) is 
given by 

1 n 

L n (Y u ..,Y n ) = (27r)-/ 2 (,;o--^n-i)- 1/2 exp(--^(y, -Ytf/v^). (7) 

i=i 

The quantity Vj = E(Yj +i — Yj + i) 2 denotes the mean-square error. It can be recursively 
calculated as well (Brockwell and Davis (1991, p. 172)). 

Next we introduce a new CUSUM chart for the variance. It is based on the idea to 
apply the likelihood ratio procedure to a Gaussian process. First, we fix n and consider 
the testing problem H Q : r > n against H± : 1 < r < n. Assuming A to be known the 
likelihood function under the null hypothesis that no change has occurred up to time n 
is equal to the joint density of (X ± , ..,X n ) in the in-control case and it is given by 

1 n 

fo(X l7 ..,X n ) = (2 7 r)-/ 2 (t; .--t; n _ 1 )- 1/2 exp(--^(X J -A > J ) 2 M-i) (8) 

i=i 

with Vj as above and X t = X]*=i a tjXj- 

If there is a change at position r e {1, ...,n} then the joint density of (Xi, ..,X n ) is 
equal to 

XX 1 
f T (X 1 , . . . , X n ) = fo(Xi, . . . , X T _i, — , . . . , — ) x 



A ' " ' A A n ~ r+1 



- 1- 

3=1 3 
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with 



7 - / V ' fOT ' ' • ' 7 V r/ V 



Xj/A for r < j < 



n 



v=l 



Let denote the indicator function of the set A and let Tj T = Yli= T a jvX v . Note that 
T jjT = for j < t. Then 



min{j-l,T-l} j-l | 

^ a^X„ + — a jvX v — Xj + (— - l)T jj7 

?J = 1 U=T 



and 



7 — 7 



X 3 " X 3 



for 1 < j < T 

X,-X 3 + ^-l){X -T hT ) for r<j<n 



Thus it follows that 



f T (X u ...,X n ) = (2 7 rr> ...^_ i r 1/2 K ^ TT exp|- 1 ^ & X ^ 



E 



(9) 



The likelihood ratio is given by 

f (Xi, ...,X n ) . . A „_ r+ i 
— — — - = mm{l, mm A 

max f T (X U ...,X n ) l<r<n 
0<r<n 



X 



exp{ -l ( j2 {Xj - Xj? - {Xj - Xj + ^ - 1){Xj - Tj - r)) ' 



\3=T 



Vj-1 



I}}- 



Hence, 



-2 log 



/o(Xx X w ) | = Q (_ (n _ T + 1)1 (A 2 } 

max f T (X 1 ,...,X n ) } 'i<r<nV 

. 0<r<n 



n 1 



2(1 - - X 3 ){X 3 - T jtT ) - (1 - ^) 2 (X, - T JiT ) 2 



Replacing A by A* > 1 the run length of the corresponding control chart is given by 



N LR (c; A*) = inf{n G IV : max{0, max (-(n - i + 1) log(A* 2 ) 

Ki<n V 



n 1 



2(1 - - X 3 ){X, - T hl ) - (1 - i^) 2 (X, - Tjj) 



}>c} 
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with c > 0. 

In the next section the control statistic is determined for several important special 
cases. Note that the above result includes all causal autoregressive moving average pro- 
cesses with normal white noise and all causal autoregressive fractionally integrated moving 
average processes with normal white noise and \d\ < 0.5. 

3.3 Examples 

In the following we will make use of the notation 



j=l J 



(10) 



with Xj as in Section 3.2. Note that in the case of an independent random sequence 
S n (A) is equal to S n (A) (see Q). 

The most popular family of time-correlated processes are autoregressive moving aver- 
age processes (ARMA). A stochastic process is called an ARMA process of order (p, q) if 
it is a solution of the stochastic difference equation 

V <i 

Y t = J2 & Y t-i + £ t + 6 i £t ~r 

i=l j=l 

Here it is assumed that {et} is an independent and normally distributed random 
process with E(e t ) = and Var(e t ) = o 2 . Moreover, it is demanded that the roots of 
1 — Yli=i $i z% are an ly m § outside the unit circle. Then the ARMA process has a unique 
stationary and causal solution. 

a) For an AR(1) process with |0i| < 1 we get that Y t = <piY t _i for t > 2 and Y\ = 0. 
Consequently we get that v t = a 2 for t > 1 and v o = 7o = o" 2 /(l — <p 2 )- Furthermore 
<h,t-i — 4>i f° r t>2, atj — for j — 1, ..,t — 2 and 2} jT = faXj-x = Xj for j > r + 1 and 
T jfT = for 1 < j < t. 

We get 



max I -(n - r + 1) log(A 2 ) + — — 

l<r<n \ V ' V ' ^ Vj-l 



'2(1 - - X 3 )(X 3 - Tj T ) - (1 - i) 2 (X, - T hT f 



max ,-(n-r + l)log(A 2 ) + (l-i,) £ ( *' * j) + 2(1 - \)^^X T - (1 - \) 2 
l<T<n \ A z Vj-i A v T -i A 

t jf=r+l J 
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;i - 7^) max (S n (A) - S T (A) + 2 — i— X T - 1— - A"( A) 

A z l<T<n \ 1 + 1/A f r _l 1 + 1/A V T -i 

;i--L) max ( S n (A) - S T (A) + - K(A) ' ArAr 



A 2 KKnl V T -i 1 + 1/A V T _ 



= (1 - h max (s n (A) - 5 T _ X (A) ^ ) . (11) 

A 2 l<r<n \ V T _x 1 + 1/A V T _i I 

One of the problems in calculating the above quantity consists in the fact that the max- 
imum has to be taken over all time points. This is usually quite time consuming and 
makes a procedures inattractive. In the present case, however, it is possible to derive a 
recursion. Let 

A+(A) = max ( S n (A) - S T ^(A) - + ^ 

Kr<n 1 U T _! 1 + 1/A f T _i 

for n > 1 and Aq(A) = then it holds for n > 1 that 

A~{A) = nmx<{ {Xn ~ ^ - K(A) - *L + 2/A 



fn-1 1 + 1/A D n _i 

max I ,S'„(A)-.SV-.l(A)-^-+ 2/A * t * t 



l<r<n-l \ W T _ 1 1 + 1/A U r - 



1 



= (X " - Xn? - K{A) + max | + ,(A) 

This representation turns out to be quite useful since it shows that the decision rule of 
the control scheme can be calculated recursively. Replacing A by A* the run length of 
this control scheme is given by 

N LR (c; A*) = inf {n G N : max{0, A+(A*)} > c] . (12) 



Putting 0! = in (11) we get the CUSUM variance chart for independent variables 



which was discussed in Section 3.2. Moreover, the representation (11) shows as well the 
relationship to residual charts. The control statistic of the CUSUM variance chart for 
independent samples applied to the residuals has a similar structure (cf. Section 4). The 
difference lies in the additional quantities based on A 2 and X T X T . 
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b) Let {Y t } be a causal AR(2) process. Then Y t = 4>iY t -i + 02^t-2 f° r t > 3, 



F 2 = 



01 



-Fi, and Yi = 0. 



1-02 

Thus we have for t > 3 that a t ^i = <p 1 , a^t-i = 02 and a t j = for 1 < j < t — 3. 
Moreover, a 2 i = 0i/(l - 02)- Furthermore, v = 7o, vi = 7 (1 - 0?/(l - 02) 2 ), v t = o 2 
for t > 2, and 



7o = cr 



1-02 



21 • 



This leads to 



Tj,r = { 



(l+0 2 )[(l-02) 2 -0?] 

Xj for j > r + 2 

/>iX r for j = t + 1 > 3 

X 2 for j = t + 1 = 2 

for 1 < j < t 



Consequently, 

n 

E- 

^— ^ 7), 



Vj-1 



2(1 - - X^Xj - T hT ) - (1 - i) 2 (X, - T hT f 



=T+1 



7J r _l 



1 + A 



12 - A - 1 

I{2,3,..,n}( T ) ~( A ii 02^r-l(^r+l " ^r+l) 



tj t v A + 1 



A + l 



and 



max 

Kr<n 



n 

. (n _ r + l)lo g (A 2 ) + ^ — 



2(1 " ^)(Xj - Xj){Xj - Tj, T ) - (1 - ^)\X 3 - T hT ) 



2A 



f 1 " a*> g?„ 1 *« (A) " *< A > " K(A) + ^ " TTa x ^ )+ 



A- 1 



10 



As in the previous example, it is possible to calculate this quantity recursively. Let 
B+(A) = max (s n (A) - S T (A) - K(A) + — (X 2 T - ^X T X T ) + 

l<r<n \ V T -i 1 + A 

then 

S+(A) = {Xn - Xn? - K(A) + max ( — - -^X n _ 1 AV 1 ) + 

V n -1 \Vn-2 1 + A 

+/ {2j3 ,.. } (n - 1) ^( ATl^X n _ 2 (X n - X n ) - |^i^ X 2_ 2) _ X(A); B +_ i(A) A . 



4 A Variance Chart based on the Sequential Proba- 
bility Ratio Test 

4.1 SPRT applied to Independent Variables 

CUSUM control charts are connected to the sequential probability ratio test (SPRT). 
Here, we derive the CUSUM procedure directly from the related SPRT. Assume that the 
variables {X t } are independent and identically normally distributed with expectation \i. 
First, we consider the simple testing problem Hq : Var(X t ) = 70 against HI : Var(X t ) = 
A* 2 7o with known A*. The SPRT says that sampling is stopped at time n if S^A*) ^ 
[A, B] with S n (A) as in If S n (A*) > B then H is rejected. Otherwise, if S n (A*) < A, 
then Hq is accepted. 

Because we are interested to detect an increase in the variance we put A = 0. Moreover, 
if S n (A*) < 0, then the chart is restarted at point zero and the procedure continues. 
Setting B = c we get the standard CUSUM chart of Section 3.1 with run length 

NsprtMc; A*) = mf{n G N : max (S n (A*) - Si (A*)) > c}. 

0<i<n 

The decision rule can be recursively calculated by using S+(A*) from Q. Note that it 
holds that N SPRTjiid (c; A*) = N LR>iid (c] A*), i.e. the likelihood ratio approach and the 
sequential probability ratio procedure lead to the same control scheme if the underlying 
process consists of independent random variables. 
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4.2 SPRT applied to Gaussian Processes 

In this section we apply the SPRT to a Gaussian process as in Section 3.2. We consider 
the simple testing problem Hq against H* (see Section 4.1) with known A*. Using ([8]), 
rt9J) and the fact that T^ T =\ — Xj it follows that 



i f fr=i(Xi, X 2 , ■ ■ ■ ,X n )\ n / A *2\ , 1 /i 1 \ (Xj — X 



" (1 - T- 

2 v A* 



with S n (A) as in (10). 



Following the procedure described in the first part the SPRT leads to a CUSUM chart 
with run length 

N SPRT (c; A*) = mi{n G N : max (S n (A*) - S^A*)) > c}. (13) 

0<i<n 

Note that the decision rule can be calculated recursively as described in ([5]). This is a 
great advantage of this scheme in comparison with the CUSUM schemes obtained by the 
likelihood approach. 

Applying the SPRT approach we get a CUSUM scheme for Gaussian processes which 
is equal to the classical CUSUM chart for independent samples obtained in Sections 3.1 
and 4.1 if the recursion is applied to the normalized residuals. Thus the chart is equal to 
the CUSUM residual chart for the variance. 

Because the normalized residuals are independent the recursive presentation shows 
that the control statistic follows a Markov process. Thus for the calculation of the ARL 
and the average delay the Markov chain approach of Brook and Evans (1972) can be 
applied. Another advantage of this scheme is based on the fact that the residuals do not 
depend on the process parameters in the in-control state and thus the control limit does 
not depend on the process parameters. Consequently this approach has some advantages 
which simplify its application. 



5 A Variance Chart based on the Shiryaev-Roberts 
Approach 

The Shiryaev-Roberts (SR) approach is based on papers of Shiryaev (1963) and Roberts 
(1966). We make use of the change point model introduced in Section 2. In Section 3 
the maximum of the likelihood ratio is taken over all possible positions of the change 
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point, i.e. r G {1, ...,n}. In the SR procedure the maximum is replaced by the sum over 
t G {l,...,n}. This procedure can be interpreted as a Bayesian procedure where r has 
a geometric prior distribution with parameter p converging to 0. Pollak (1985) proved 
for independent variables that the SR-rule to be asymptotically Bayes risk efficient as 
p — > 0. The SR approach for independent variables has been recently discussed by several 
authors, e.g., Moustakides et al. (2009) and Pollak and Tartakovsky (2009). 



5.1 Shiryaev- Roberts Approach for Independent Variables 

First it is assumed that variables {Y t } are independent and normally distributed with 
mean p and variance 70 . Let / denote the density of a univariate normal distribution 
with mean p and variance 70 and let g be the density of a univariate normal distribution 
with mean p and variance A 2 7 . Using the notation of Section 3 and Lj(x) = g(x)/ f(x) 
the SR statistic is given by 

Rn(A) = J2 {fx''' )X x\ = J2fl L i = ( 1 + Rn-^))K (14) 

i=l n > i= i j =i 

for n > 1 and Ro(A) = 0. For normal variables we get that 

1 f 1 1 w„ ,2 



= (l + it^AVj-expj— (1-— )(X n -pYj. (15) 

However, in practice, we do not know the size of the change A. As described above it is 
replaced by a known quantity A* > 1. This leads to R n (A*). 

For deriving the control statistic it was assumed that the target process is independent. 
Following Section 3.1 we apply this statistic for time series as well. The independent 
variables are replaced by the time series. The run length of the chart is given by 

N SR , td (c; A*) = inf{n G N : R n (A*) > c}. (16) 
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5.2 Shiryaev- Roberts Approach for Gaussian Processes 

Now we apply this approach to a Gaussian process {Y t } fulfilling the assumptions of 
Section 3.2. Using (j8j) and ^ we receive that 

RJA) = V fA^ll • • • 

~[ fo(Xi, . . . ,X n ) 

_ \ 1 A {Xj - x 3 Hj- 1){X 3 - T 3 ,)f 1 " (X 3 - x 3 f 



A n-i+l r 2 ^ V, ! 2 ^ U, ! 

i=l ^ j'=i J j=i J 

- 1 f ™ 1 / 1 11 \ 

= E 6XP E — U 1 " a )(X ^' - - T - } - 2 (1 " A )2(X ' - T - )2 

i=i I j=i 3 1 \ ' 

Replacing A by A* the run length of the SR chart is given by 

N SR (c; A*) = inf{n e N : i? n (A*) > c}. (17) 

5.3 Example 

Suppose that {Y t } is a causal AR(1) process. Following Example a) of Section 3.3 we get 
that 



2 XiXi 



i=l (. \ j=i J 

It is possible to determine the SR statistic recursively 

MA) = ( Rn-i(A) + exp {(1 - 1 )(-J_^ - JL) 1U exp {1(1 1 ^ (X " - X " 



A 2/V l + A Wn _! 2v n H ] ) A 1 \2 K A 2/ v n _! 

for n > 1 and i?o(A) = 0. As above the unknown magnitude of the change A is replaced 
by a reference value A*. This leads to R n (A*). 

6 A Variance Chart based on the Generalized Like- 
lihood Ratio Approach 

The disadvantage of the procedures considered in Sections 3 to 5 consists in the fact that 
for the derivation of the procedures the magnitude of the change has to be known. In 
practice, however, in many cases no information is available about the size of a possible 
shift. In that case other procedures must be favored. Such approaches are discussed 
in Sections 6 to 8. In this section we consider the generalized likelihood ratio (GLR) 
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approach which can be considered as an extension of the likelihood ratio method since 
the size of the change is assumed to be unknown. GLR charts have not received much 
attention in SPC up to now (cf. Reynolds and Lou (2012)) and have been mostly discussed 
in literature on change-point analysis (e.g., Lai (2001)). 

6.1 GLR applied to Independent Variables 

Assume that {Y t } is an independent normally distributed random sequence with mean ji 
and variance a 2 . We make use of the change point model introduced in Section 2. Let 
fr,A denote the out-of-control density of (X 1: ..,X n )' then we have to calculate 

max suplog(/ Ti Apfi, ...,X n )). 



l<T<n 



A>1 



We get with f n = £" =1 pQ - /j) 2 / 1o that 



n, _ . n-T + 1, . . 9 . 1 



2 



\og(f T , A (Xi,...,X n )) = --log(27r 7 o) log(A 2 ) - - [T T _! + (T n -T T ^)/A j 

Let A 2 Tjn = (f n - f T _i)/(n - r + 1) and let A 2 T , a = max{l, A 2 n }. It holds that 
suplog(/ T)A (X 1 , ...,X n )) = log(f rA (X u ...,X n )) 

A>1 

since the derivative of log(f TA (Xi, X n )) with respect to A 2 is positive if A 2 < A 2 r 
and it is negative for A 2 > A 2 n . It holds that 

log(f rAT jX 1: ...,X n )) 

-f log(27r 7 o)-|[T r _ 1 + n-r + l]-^±ilog(A 2 n ) if A 2 n > 1 
-flog(27r 7o )-f if A 2 <1 • 



Thus 



/ .a(-V;. X n ) \ n - T + 1 2 2 



sup log """T = 1 — ^[A.,n - 1 - log(A 2 ,J]. 

A>1 V M A 1> •••) A nJ / ^ 



Because x — 1 — log(x) > for a; > 1 we obtain that 

2 max sup log / ^(X,, X ) \ = _ ^ + _ x _ 

l<r<n+l A>1 y / (Ai, A n J / l<r<n 

Consequently the stopping rule of the GLR test is given by 

N GLR , iid {c) = inf{n e N : max {(n - % + 1) [A 2 n - 1 - log(A 2 J] } > c}. 

K't<n 
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6.2 GLR applied to Gaussian Processes 

Suppose that {Y t } is a Gaussian process with mean zero and covariance function k(i,j) = 
Cov(Xi,Xj) fulfilling the conditions of Section 3.2. Using T n = J27=i(-^i ~~ ^i) 2 / v i-i an d 
Q it follows that 

log^Apf!,...,^)) = const log(A 2 ) 

"2 I Tt ~ 1 + 4- ^ 

Now we will maximize the function log / t ,a(^i, • • • ,X n ) with respect to A. For that 
reason we calculate the derivative 

^iog(/ T ,A(x l5 ...,x n )) ^ + ^2^ — (** 

7l - r + 1 + J_ (o _c ^ + J_5 



where 



2 



q _ ( X j ~ Xj)(Xj ~ T jtT ) » _ (jfj ~ ?j,T, 

Let 

S n ,r — S n ,T + \J (S njT — S n>T ) 2 + 4(n — r + l)^^ 
Ar ' n = 2(n-r + l) ' 

Because the derivative of log(/ Tj A(^i, ...,X n )) with respect to A is positive if < A < A T)fi 
and it is negative else, it holds that 



suplog(/ TjA (X 1 ,...,X n )) = log(/ T)A * n CXi,...,-Xn)) 

A>1 

with A* n = max{l, A T)n }. Thus 

, ( fr,A(Xl, ...,X n )\ 

sup log ~TJy YT 

A>1 V M A 1> A n J / 

= _(„ _ r + i) log(A;j - - l)(2^, r + (-L - l)S n , T ). (19) 

Consequently the stopping rule for this chart is given by the following 
Nglr(c) = inf \ n G N : max { - (n - i + 1) log(A*J 

I l<i<n K 
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6.3 Example 

Let {Y t } be a causal AR(1) process. Using the results of Section 3.3a) we get that 

S nT = y & - - Tj > } =T n -T T+ {Xt - * t)Xt , (20) 

~ Vj-l ' Ur-l 

S n , T = j2 iXj ~ Tj ' T)2 =T n -T T + ^. (21) 
Thus S n>T = S UtT — x v tX ^ ■ Consequently 



S n ,r — S n:T + y (Sn !T — Sn :T ) 2 + — T + 1)5^,- 

2(n - r + 1) 

X T X T | , I( X T X T \ 2 i 1L_ T x 1VT _T _l_ 



= ^ + V(^) 2 + 4 ^-^+ 1 )( T »-^ + ^ : 

2(n - r + 1) 

and the stopping rule is given by 

N G lr{c) =inf(nGN: max {-(n - i + 1) log(A* n ) - 1 ' 



2 V A*„ 

1 +l)(T n -T l + ^)-2^)>cl. 



i,n V i-1 V i~l 

7 A Variance Chart based on the Generalized SPRT 

Let {Y t } be a Gaussian process as assumed in Section 3.2. Following the SPRT approach 



sketched in Section 4 and using (19) we get 



sup log 



/r=l,ApG.) •••) X n ) 



A>1 V fo(X\, ...,X n ) 



-nlog(Ay - - l)(2S n> i + (^r - 



17 



Since S U: i = S Uj i = T n we get that Ai jn = \jT n /n and thus 

, f /t=i,a(-Xi, X n ) 
sup log ~~ TTy yY~ 

A>1 V J0{^1, •••) -^n) 

-n [log(T„/n) - T n /n + 1] /2 if T n /n > 1 

n -f t / / i = ^( T -/ n ) 

it T n /n < 1 

with 

fc n (a?) = nh(x) =n(x-l- log(x))/2. (22) 

Note that h n (T n /n) > 0. Following Section 4 a control chart is obtained by applying this 
approach sequentially. The run length of this scheme is 

Ngsprt(c) = inf{n e N : max (h n (T n /n)) - h^Ti/i)) > c}. 

0<i<n 



Assuming r = 1 in (Fn) it holds that P T= \^(T n < x) = %^(x/A 2 ) and thus P T =i : A(T n /n < 
1) = Xn( n /A 2 ) — > if A — > oo. This means that the probability that h n (T n /n) is positive 
is increasing with A. 

8 Variance Charts based on the Generalized Shiryaev- 
Roberts Approach 

8.1 GSR applied to Independent Variables 

Assume that the variables {Y t } are independent and normally distributed with mean fi 



and variance 70. Following (14) and (15) it is necessary to determine the maximize of 
R n (A) over A. However, it is not possible to get an explicit expression for the value of 
A which maximizes this quantity since the derivation of i? n (A) with respect to A leads 
to an exponential sum which is difficult to handle. For that reason we choose another 
procedure. Instead of -R n (A) we consider 



D*/A\ V^l f fi,A{Xi, ...,X n )\ 



Because the logarithm is a strictly increasing continuous function this means that instead 
of the arithmetic mean the geometric mean is maximized. 
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In the case of independent variables we get that 

D* (A\ 1 I TT /*,a(-^1) Xn) \ 

fl„,«,(A) = l°g(TI WXl „„ iA y J 

= £ (-(" - * + 1) 'os(A) + (i - ^)(r„ - ■ 

Determining the derivative of R* niid (A) we see that the maximum of R* niid (A) is attained 
at 

n(n + 1) 

with 

— 7o 

i=l i=l ,u 

Because i?* iid (A) is a concave function the maximum over A > 1 is attained at A n ^d = 
max{l, A ntUd }. For A niid < 1 it follows that supi?* iid (A) = 0, else 

A>1 

2 sup K,^(A) = ~ 1 9 ^ log(A^ ) + (1 - l/A^ iid )C/ n 

A>1 ^ 

n(n + l) n(n+l) / , . fn(n+l) 
= U n I \og{U n ) - log 



n(n + 1) 



(A 



nAid! 



with h n as in (22). The run length of the control chart is given by 

N G SR,ud(c) = inf{n E N : /i n(n+ i)(A^ iid ) > c}. 

8.2 GSR for Gaussian Processes 

Now let {Yj} be a Gaussian process as assumed in Section 3.2. As in the previous section 
we make use of the geometric mean i?* (A). In the present case we get that 

2i?;(A) = ^L(n-A; + l)log(A 2 ) + 2(l-i)5 n , fc -(l-i) 2 5 n , fc ') 
k=i ^ ' 
n(n + 1) /a9x . 1 . ■ , 1 ,n 



log(A 2 ) + 2(1 - -)U n - (1 - -)^ n 



2 " v 7 A' " v A 
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with 



fc=l 



k=l 



and S n ,k and S^ft as in (18). R* n (A) is a concave function for A > and its maximum for 



A > is attained at position 



Un - U n + sJ{U n -U n y + 2n{n + l)U n 



n(n + 1) 



Thus it holds that 



2supK(A) 

A>1 



n(n + I") 



log(A') + 2(1 - j-ft 



A, 



i U n = g n (U n ,U n ) 



where A n = max{l, A n }. The run length of this scheme is 



Ngsr(c) = inf{n E N : g n (U n , U n ) > c}. 



9 Comparison Study 

In the above sections several new control schemes for detecting a change in the variance of 
a time series were introduced. Because no optimality results are known our aim is to give a 
practitioner a clear recommendation which chart should be applied in a specific situation. 
In order to compare control charts at all, we have to calibrate them. The control limits 
of all charts are chosen such that their in-control ARLs are the same. Here we fix the 
in-control ARL to be 500. Using the corresponding control limits the out-of-control ARLs 
of all charts are compared with each other using the out-of-control ARL E T= i^(N(c)) as 
well as the average delay E Tjl \(N(c) — r + l\N(c) > r). These performance criteria are 
mostly applied in literature. 

In the present case we do not have explicit formulas for the calculation of the per- 
formance criteria. Such formulas are not available for control charts for time-dependent 
processes. Here we make use of an extensive simulation study. In each case the average 
run length and the average delay were determined within a simulation study based on 10 6 
repetitions. The only exception are the GLR charts where no recursive presentation was 
given and the calculation of the average run length and the average delay is more time 
consuming. In that case we used 10 5 repetitions. 

Note that some control charts depend on a reference value A* > 1. In our study A* 
is taking values within the set {1.10, 1.20, 1.30, 1.40, 1.50, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00}. 

In our comparison study the target process is an AR(1) process with standard nor- 



20 



mally distributed white noise. The coefficients of the process are taking values within 
{-0.9, -0.8, -0.7, . . . , 0.7, 0.8, 0.9}. 

In Tables 1 to 3 the out-of-control ARLs of the considered charts are given. Because 
the ARL turns out to be symmetric with respect to the coefficient (j>\ of the AR process 
we only show the results for non-negative values of 0i. In each row and each column the 
ARLs of nine control charts are given, above the variance chart for iid variables (Section 



3.1, cf. (§)), followed by the LR chart (Section 3.2, cf. (12)), the SPRT chart (Section 4.2 



cf. (13)), the Shiryaev-Roberts chart for iid variables applied to time series (Section 5.1, 



cf. (16)), the Shiryaev-Roberts chart for Gaussian processes (Section 5.2, cf. (17)), the 
GLR chart of Section 6.2, the GSPRT chart of Section 7, and the GSR chart of Section 
8.2. The first five charts depend on a reference value. For these charts the smallest out- 
of-control ARL over all A* is listed. In parenthesis the value of A* is given where the 
minimum is attained. The other four charts are generalized schemes and do not depend 
on a reference value. In Tables 1 to 3 the ARLs of all charts are written in bold which 
for a fixed value of A deviate from the smallest out-of-control ARL by only 2%. 

[ Tables 1 to 3 about here. ] 

The results of the comparison study are very interesting. First, it can be seen that the 
results for the charts based on the independence assumption, i.e. the schemes of Section 
3.1 and 5.1, are getting worse if the correlation structure of the target process increases. 
Since the other schemes behave much better they should not be applied. Second, the 
minimum out-of-control ARL of the LR chart and the SPRT chart is always attained 
if A* is equal to the true value of the change A. For the SR chart a slightly different 
behavior is observed. Here the best value is greater or equal to A. If the optimal values 
for the LR and the SPRT chart are taken they provide smaller ARLs than the SR chart. 
Among the five schemes with a reference value the LR and the SPRT chart behave the 
best. Their results are very similar. 

The analysis of the generalized charts shows that again the chart based on the indepen- 
dence assumption (here: GSRiid, Section 8.1) behaves bad. The results for the GSPRT 
chart are also not very good for small changes but it is the best generalized scheme for 
A > 1.75. The best results for small changes (A < 1.3) are obtained for the GSR chart. 
The chart even behaves better than all charts using a reference parameter. This is very 
remarkable. However, for A > 1.4 the LR chart and the SPRT chart dominate this 
scheme. 

In Figure 1 we discuss the sensitivity of the LR, the SPRT, and the SR chart with 
respect to the choice of the reference value. In practice we usually do not know the 
magnitude of the expected change in advance. In the figure the results of the GLR, the 
GSPRT, and the GSR chart are given as well. We focus on two changes, A = 1.3 and 
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A = 2.0 and on two values of the coefficient of the AR(1) process, 0i = 0.4 and <pi = 0.8. 
The figure shows the dominance of the GSR chart for A = 1.3. For a small correlation 
structure, here 0i = 0.4, the GLR chart and the SPRT chart show to be better than the 
GLR chart if A* < 2.0. This means that if the deviation from the optimal A = 1.3 is not 
too large then these charts have a smaller ARL than the GLR chart but if the deviation 
is large (A* > 2.25) then the GLR scheme must be preferred. A similar behavior can be 
observed for (pi = 0.8. However, if the change is larger, here A = 2.0, then the LR and the 
SPRT chart are always better than the GLR chart. Even if the reference value is chosen 
completely different than the true value of the change the schemes are better. Only the 
GSPRT chart turns out to be better than the LR and the SPRT scheme if the deviation 
from the optimal value is sufficiently large. It is interesting as well that it seems to be 
better to overestimate the value of the change than to underestimate it. 

[ Figure 1 about here. ] 

Up to now the charts were compared using the average run length. For this perfor- 
mance measure it is assumed that the change already arises at the beginning, i.e. r = 1. 
This is of course a restriction. The average delay is a more general criteria because the 
change may arise at any position. In Table 4 it is assumed that the change arises up to 
the 50th observation, i.e. 1 < r < 50. Moreover, we focus on the changes A = 1.3 and 
A = 2.0. The values of the LR, the SPRT, and the SR chart refer to the optimal choice of 
the reference value. The table shows that except for the GSPRT chart the worst average 
delay is attained at r = 1, i.e. it is equal to the ARL. As r increases the average delay 
is decreasing and it does not change a lot for r > 20. The GSPRT chart is an exception. 
The worst average delay is attained at r = 50 and its minimum ARL is observed for a 
small value of r. The table illustrates that the GSR chart behaves quit well for small 
changes. The limit of the average delay and the worst case average delay are the smallest 
ones for this chart. However, for large changes its behavior is more complicate. While 
its ARL is the largest one it turns out to be more effective if the change arises at a later 
time point. 

[ Table 4 about here. ] 

10 Summary 

In this paper several new control charts for the detection of a change in a Gaussian process 
are introduced. The charts are derived by making use of the likelihood ratio approach, the 
sequential probability ratio test, and the Shiryaev- Roberts approach. For the derivation 
of the charts it is assumed that the size of the change is known. The obtained charts 
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depend on a reference value which has to be chosen suitably. We consider the case of an 
unknown size of the change as well. This attempt leads to generalized control charts. 

In an extensive simulation study we compare the introduced control charts with each 
other. The target process is assumed to be an AR(1) process. Using the ARL as a 
performance criterion it turns out that the GSR chart behaves the best for small changes 
(A < 1.3). For detecting medium and large changes (A > 1.3) the LR and the SPRT 
chart turn out to be better. They depend on an additional reference value. The minimum 
out-of-control ARL is obtained for both schemes if the reference value is chosen equal to 
the true value of the change. It turns out that for larger changes the LR and the SPRT 
chart are still better than the best generalized chart if the true change is not dramatically 
underestimated. 

If we analyze the control charts using the average delay it can be seen that except the 
GSPRT chart the worst average delay is always attained at r = 1, i.e. it is equal to the 
ARL. The GSR chart provides the best results for a small change. For medium and large 
change the LR and the SPRT chart must be favored if the change is expected to be at 
the beginning of the monitoring process. However, if it appears at a later time point the 
GSR turns out to have the smallest delay. 
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Figure 1: Out-of-control ARLs of several CUSUM charts as a function of the reference 
value A* for an in-control ARL of 500 
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Table 4: Average run length (above), worst average delay for 1 < r < 50 (middle), and 
the value of the delay at position r = 50 (below) for the LR, the SPRT, and the SR chart 
for optimal reference parameter, and the GLR, the GSPRT chart, and the GSR chart 



(01 


= 0.4, 


in-control ARL = 


= 500) 










LR 


SPRT 


SR 


GLR 


GSPRT 


GSR 


A 


= 1.3 


32.52 


32.59 


35.09 


39.40 


49.18 


32.68 






32.52 


32.59 


35.09 


39.40 


73.46 


32.68 






29.85 


29.85 


30.70 


33.43 


73.46 


21.91 


A 


= 2.0 


7.52 


6.79 


7.14 


9.42 


8.38 


11.41 






7.52 


6.79 


7.14 


9.42 


14.22 


11.41 






6.69 


6.41 


6.48 


8.16 


14.22 


6.18 
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